C*************************************************************** C Subroutine LSEM3 by Stephen Kirkup C*************************************************************** C C Version 2 (July 2015) C New GLS routine implemented and simplified integral equation C formulation C Stephen Kirkup, School of Engineering C University of Central Lancashire - www.uclan.ac.uk C smkirkup@uclan.ac.uk C http://www.researchgate.net/profile/Stephen_Kirkup C C This open source code can be found at C www.boundary-element-method.com/fortran/LSEM3.FOR C C Issued under the GNU General Public License 2007, see gpl.txt C C Part of the the author's open source BEM packages. C All codes and manuals can be downloaded from C www.boundary-element-method.com C C C*************************************************************** C This subroutine computes the solution of Laplace's equation C __ 2 C \/ {\phi} = 0 C C in the region exterior to a set of thin surfaces. C C C How to use the subroutine C ------------------------- C C The following diagram shows how the subroutine is to be used. C C .................................... C : : C : : C ---------------------- : -------------------------- : C | | : | | : C | MAIN PROGRAM |------->-----| LSEM3 | : C | (e.g. lsem3_t.for) | : | | : C | | : -------------------------- : C ---------------------- : | : C : > : C : | : C : ------------------------ : C Package ---------------->: | subordinate routines | : C : ------------------------ : C : : C : (this file) : C :..................................: C / | | C |_ > > C / | | C ................ ................ ................ C : : : -------- : : --------- : C : (geom3d.for) :---<---: | L3LC | : : | GLS | : C : : : -------- : : --------- : C :..............: : -------------: : -------------: C : |subordinate|: : |subordinate|: C : | routines |: : | routines |: C : -------------: : -------------: C : : : : C : (l3lc.for) : : (gls.for) : C :..............: :..............: C C C The contents of the main program must be linked to LSEM3.FOR, L3LC.FOR C GLS.FOR and GEOM3D.FOR. C C Method of solution C ------------------ C C In the main program, the shell(s) must be described as a set of C panels. The panels are defined by three indices (integers) which C label a node or vertex on the shell. The data structure VERTEX C lists and enumerates the coordinates of the vertices, the data C structure PIELV defines each panel by indicating the labels for C the three nodes that are its vertices and hence enumerates the C panels. C C Format and parameters of the subroutine C --------------------------------------- C C The subroutine has the form C C SUBROUTINE LSEM3(MAXNV,NV,VERTEX,MAXNPI,NPI,PIELV, C * MAXNPE,NPE,PEXT, C * HA,HB,HF,HAA,HBB,HFF, C * HIPHI,HIVEL,PEIPHI, C * LSOL,LVALID,EGEOM, C * PHIDIF,PHIAV,VELDIF,VELAV,PEIPHI, C * AMAT,BMAT,L_EH,M_EH, C * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3) C The parameters to the subroutine C ================================ C Geometry of the shell {\Pi} (input) C integer MAXNV: The limit on the number of vertices of the triangles C that defines (approximates) {\Pi}. MAXNV>=3. C integer NV: The number of vertices on {\Pi}. 3<=NV<=MAXNV. C real VERTEX(MAXNV,3): The coordinates of the vertices. VERTEX(i,1), C VERTEX(i,2), VERTEX(i,3) are the x,y,z coordinates of the i-th C vertex. C integer MAXNPI: The limit on the number of panels describing {\Pi}. C MAXNPI>=1. C integer NPI: The number of panels describing {\Pi}. 1<=NPI<=MAXNPI. C integer PIELV(MAXNPI,3): The indices of the three vertices defining C each panel. The i-th panel has vertices C (VERTEX(PIELV(i,1),1),VERTEX(PIELV(i,1),2)),VERTEX(PIELV(i,1),3)), C (VERTEX(PIELV(i,2),1),VERTEX(PIELV(i,2),2)),VERTEX(PIELV(i,2),3)) and C (VERTEX(PIELV(i,3),1),VERTEX(PIELV(i,3),2)),VERTEX(PIELV(i,3),3)). C Exterior points at which the solution is to be observed (input) C integer MAXNPE: Limit on the number of points exterior to the C shell. MAXNPE>=1. C integer NPE: The number of exterior points. 0<=NPE<=MAXNPE. C real PEXT(MAXNPE,3). The coordinates of the exterior point. C PEXT(i,1),PEXT(i,2),PEXT(i,3) are the x,y,z coordinates of the i-th C point. C Shell Conditions C a(p){\delta}(p)+b(p){\nu}(p)=f(p) C real HA(MAXNPI). Function a at the central points C real HB(MAXNPI). Function b at the central points C real HF(MAXNPI). Function f at the central points C A(p){\Phi}(p)+B(p)V(p)=F(p) C real HAA(MAXNPI). Function A at the central points C real HBB(MAXNPI). Function B at the central points C real HFF(MAXNPI). Function F at the central points C Incident field C real HIPHI(MAXNPI). The incident potential {\phi} at the central C points C real HIVEL(MAXNPI). The derivative of the incident potential {\phi} C at the central points C real PEIPHI(MAXNPE). The incident potential {\phi} at the exterior C points C Validation and control parameters (input) C logical LSOL: A switch to control whether the particular solution is C required. C logical LVALID: A switch to enable the choice of checking of C subroutine parameters. C real EGEOM: The maximum absolute error in the parameters that C describe the geometry. C Solution (output) C real SPHI(MAXNPI): The potential ({\phi}) at the centres of the C shell panels. C real SVEL(MAXNPI): The (v or d{\phi}/dn where n is the outward normal C to the shell) at the centres of the shell panels. C real PEIPHI(MAXNPE): The potential ({\phi}) at the exterior points. C Working space C Discrete Operators (AMAT and BMAT changed in GLS) C real AMAT(2*MAXNPI,2*MAXNPI) C real BMAT(2*MAXNPI,2*MAXNPI) C real L_EH(MAXNPE,MAXNPI) C real M_EH(MAXNPE,MAXNPI) C Further information from system solution GLS C integer PERM(MAXNPI) C logical XORY(MAXNPI) C real C(MAXNPI) C real HALPHA(2*MAXNPI) C real HBETA(2*MAXNPI) C real HFFF(2*MAXNPI) C real X(2*MAXNPI) C real Y(2*MAXNPI) C Work Space C real WKSPC1(MAXNPI) C real WKSPC2(MAXNPI) C real WKSPC3(MAXNPI) C Notes on the geometric parameters C --------------------------------- C (1) The indices of the nodes listed in PIELV must be such that they C are ordered counter-clockwise around the shell, when viewed C from above the shell on the + side. C Notes on the exterior points C ---------------------------- C (1) The points in PEXT should not lie on the shell, as defined C by the parameters VERTEX and PIELV. Any point lying outside the C shell. C The subroutine SUBROUTINE LSEM3(MAXNV,NV,VERTEX,MAXNPI,NPI,PIELV, * MAXNPE,NPE,PEXT, * HA,HB,HF,HAA,HBB,HFF, * HIPHI,HIVEL,PEIPHI, * LSOL,LVALID,EGEOM, * PHIDIF,PHIAV,VELDIF,VELAV,PEPHI, * AMAT,BMAT,L_EH,M_EH, * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3) PARAMETER (MAXNQ=100) C shell geometry C Limit on the number of vertices on {\Pi} INTEGER MAXNV C The number of vertices on {\Pi} INTEGER NV C The coordinates of the vertices on {\Pi} REAL*8 VERTEX(MAXNV,3) C Limit on the number of panels describing {\Pi} INTEGER MAXNPI C The number of panels describing {\Pi} INTEGER NPI C The indices of the vertices describing each panel INTEGER PIELV(MAXNPI,3) C Exterior points at which the solution is to be observed C Limit on the number of points exterior to the shell where C solution is sought INTEGER MAXNPE C The number of exterior points INTEGER NPE C Coordinates of the exterior points REAL*8 PEXT(MAXNPE,3) C Shell Conditions C a(p){\delta}(p)+b(p){\nu}(p)=f(p) C function a at the central points REAL*8 HA(MAXNPI) C function b at the central points REAL*8 HB(MAXNPI) C function f at the central points REAL*8 HF(MAXNPI) C A(p){\Phi}(p)+B(p)V(p)=F(p) C function A at the central points REAL*8 HAA(MAXNPI) C function B at the central points REAL*8 HBB(MAXNPI) C function F at the central points REAL*8 HFF(MAXNPI) C Incident field C The incident potential {\phi} at the central points REAL*8 HIPHI(MAXNPI) C The derivative of the incident potential {\phi} at the C central points REAL*8 HIVEL(MAXNPI) C The incident potential {\phi} at the exterior points REAL*8 PEIPHI(MAXNPE) C Validation and control parameters LOGICAL LSOL LOGICAL LVALID REAL*8 EGEOM C Solution REAL*8 PHIDIF(MAXNPI) REAL*8 PHIAV(MAXNPI) REAL*8 VELDIF(MAXNPI) REAL*8 VELAV(MAXNPI) REAL*8 PEPHI(MAXNPE) C Working space REAL*8 AMAT(2*MAXNPI,2*MAXNPI) REAL*8 BMAT(2*MAXNPI,2*MAXNPI) REAL*8 L_EH(MAXNPE,MAXNPI) REAL*8 M_EH(MAXNPE,MAXNPI) C Further information from system solution GLS INTEGER*4 PERM(2*MAXNPI) LOGICAL XORY(2*MAXNPI) REAL*8 C(2*MAXNPI) REAL*8 HALPHA(2*MAXNPI) REAL*8 HBETA(2*MAXNPI) REAL*8 HFFF(2*MAXNPI) REAL*8 PHIINF(2*MAXNPI) REAL*8 VELINF(2*MAXNPI) C Work Space REAL*8 WKSPC1(2*MAXNPI) REAL*8 WKSPC2(2*MAXNPI) REAL*8 WKSPC3(2*MAXNPI) REAL*8 DIST3,AREA REAL*8 DIAM C Constants C Real scalars: 0, 1, 2, half, pi REAL*8 ZERO,ONE,TWO,THREE,THIRD C Geometrical description of the shell C panels counter INTEGER IPI,JPI C The points exterior to the shell where the solution is sought INTEGER IPE C Parameters for L3LC REAL*8 P(3),PA(3),PB(3),PC(3),QA(3),QB(3),QC(3),VECP(3) LOGICAL LPONEL C Quadrature rule information C [Note that in this program two quadrature rules are used: one for C the case when the point P lies on the panel (LPONEL=.TRUE.) and C one for the case when P does not lie on the panel. In general, C it is more efficient to define a larger set of quadrature rules C so that a particular rule can be selected for any given point P C and panel QA-QB-QC. For example using more quadrature points when C the panel is large, less when the panel is small, more when C the panel is close to P, less when it is far from P.] C Quadrature rule used when LPONEL=.TRUE. C Number of quadrature points INTEGER NQON C x-Abscissae of the actual quadrature rule REAL*8 XQON(MAXNQ) C y-Abscissae of the actual quadrature rule REAL*8 YQON(MAXNQ) C Weights of the actual quadrature rule REAL*8 WQON(MAXNQ) C Quadrature rule used when LPONEL=.FALSE. C Number of quadrature points INTEGER NQOFF C x-Abscissae of the actual quadrature rule REAL*8 XQOFF(MAXNQ) C y-Abscissae of the actual quadrature rule REAL*8 YQOFF(MAXNQ) C Weights of the actual quadrature rule REAL*8 WQOFF(MAXNQ) C Quadrature rule parameters for L3LC C Actual number of quadrature points INTEGER NQ C Abscissae of the actual quadrature rule REAL*8 XQ(MAXNQ) C Abscissae of the actual quadrature rule REAL*8 YQ(MAXNQ) C Weights of the actual quadrature rule REAL*8 WQ(MAXNQ) C Counter through the quadrature points INTEGER IQ C Validation and control parameters for subroutine L3LC LOGICAL LVAL REAL*8 EQRULE LOGICAL LFAIL1 LOGICAL LL LOGICAL LM LOGICAL LMT LOGICAL LN C Parameters for subroutine L3LC. REAL*8 DISL REAL*8 DISM REAL*8 DISMT REAL*8 DISN C Other variables C Failure flag LOGICAL LFAIL C Accumulation of solution {\phi} REAL*8 SUMPHI C Error flag LOGICAL LERROR C INITIALISATION C -------------- C Set constants ZERO=0.0D0 ONE=1.0D0 TWO=2.0D0 THREE=3.0D0 THIRD=ONE/THREE C Set up validation and control parameters C Switch off the validation of L3LC LVAL=.FALSE. C Set EQRULE EQRULE=1.0D-6 C Set EGEOM EGEOM=1.0D-6 C Set up the quadrature rule(s). C Set up quadrature rule for the case when P is not on the panel. C Set up 7 point Gauss-Legendre rules CALL GLT7(MAXNQ,NQOFF,WQOFF,XQOFF,YQOFF) C Set up quadrature rule for the case when P is on the panel. C Set up quadrature rule data. If LPONEL is false then use the standard C Gaussian quadrature rule above. If LPONEL is true then then a C quadrature rule with 3 times as many points is used, this is made C up from three standard quadrature rules with the quadrature points C translated to the three triangles that each have the cetroid and two C of the original vertices as its vertices. NQON=3*NQOFF DO 330 IQ=1,NQOFF XQON(IQ)=XQOFF(IQ)*THIRD+YQOFF(IQ) YQON(IQ)=XQOFF(IQ)*THIRD WQON(IQ)=WQOFF(IQ)/THREE XQON(IQ+NQOFF)=XQOFF(IQ)*THIRD YQON(IQ+NQOFF)=XQOFF(IQ)*THIRD+YQOFF(IQ) WQON(IQ+NQOFF)=WQOFF(IQ)/THREE XQON(IQ+2*NQOFF)=THIRD*(ONE+TWO*XQOFF(IQ)-YQOFF(IQ)) YQON(IQ+2*NQOFF)=THIRD*(ONE-XQOFF(IQ)+TWO*YQOFF(IQ)) WQON(IQ+2*NQOFF)=WQOFF(IQ)/THREE 330 CONTINUE C Validation C ========== C Validation of parameters of LEBEM3 C --------------------------------- IF (LVALID) THEN C Validation of main paramters LERROR=.FALSE. IF (MAXNV.LT.3) THEN WRITE(*,*) 'MAXNV = ',MAXNV WRITE(*,*) 'ERROR(LSEM3) - must have MAXNV>=4' LERROR=.TRUE. END IF IF (NV.LT.3.OR.NV.GT.MAXNV) THEN WRITE(*,*) 'NV = ',NV WRITE(*,*) 'ERROR(LSEM3) - must have 3<=NV<=MAXNV' LERROR=.TRUE. END IF IF (MAXNPI.LT.3) THEN WRITE(*,*) 'MAXNPI = ',MAXNPI WRITE(*,*) 'ERROR(LSEM3) - must have MAXNPI>=3' LERROR=.TRUE. END IF IF (NPI.LT.3.OR.NPI.GT.MAXNPI) THEN WRITE(*,*) 'NPI = ',NPI WRITE(*,*) 'ERROR(LSEM3) - must have 3<=NPI<=MAXNPI' LERROR=.TRUE. END IF IF (MAXNPE.LT.1) THEN WRITE(*,*) 'MAXNPE = ',MAXNPE WRITE(*,*) 'ERROR(LSEM3) - must have MAXNPE>=1' LERROR=.TRUE. END IF IF (NPE.LT.0.OR.NPE.GT.MAXNPE) THEN WRITE(*,*) 'NPE = ',NPE WRITE(*,*) 'ERROR(LSEM3) - must have 3<=NPE<=MAXNPE' LERROR=.TRUE. END IF IF (EGEOM.LE.ZERO) THEN WRITE(*,*) 'NPE = ',NPE WRITE(*,*) 'ERROR(LSEM3) - EGEOM must be positive' LERROR=.TRUE. END IF IF (LERROR) THEN LFAIL=.TRUE. WRITE(*,*) WRITE(*,*) 'Error(s) found in the main parameters of LSEM3' WRITE(*,*) 'Execution terminated' STOP END IF END IF C Find the diameter DIAM of the shell DIAM=0.0 DO 100 IV=1,NV-1 PA(1)=VERTEX(IV,1) PA(2)=VERTEX(IV,2) PA(3)=VERTEX(IV,3) DO 110 JV=IV+1,NV PB(1)=VERTEX(JV,1) PB(2)=VERTEX(JV,2) PB(3)=VERTEX(JV,3) DIAM=MAX(DIAM,DIST3(PA,PB)) 110 CONTINUE 100 CONTINUE IF (LVALID) THEN LERROR=.FALSE. C Check that the shell defined by PIELV is complete and closed DO 120 ISE=1,NPI 120 CONTINUE C Check that EGEOM is not too large IF (EGEOM.GT.DIAM/100.0D0) THEN WRITE(*,*) 'EGEOM = ',EGEOM WRITE(*,*) 'ERROR(LSEM3) - EGEOM is set too large' LERROR=.TRUE. END IF IF (LERROR) THEN LFAIL=.TRUE. WRITE(*,*) WRITE(*,*) 'Error in shell geometry or EGEOM' WRITE(*,*) 'Execution terminated' END IF END IF IF (LVALID) THEN C Check that the vertices are distinct with respect to EGEOM LERROR=.FALSE. DO 130 IV=1,NV-1 PA(1)=VERTEX(IV,1) PA(2)=VERTEX(IV,2) PA(3)=VERTEX(IV,3) DO 140 JV=IV+1,NV PB(1)=VERTEX(JV,1) PB(2)=VERTEX(JV,2) PB(3)=VERTEX(JV,3) IF (ABS(PA(1)-PB(1)).LT.EGEOM) THEN IF (ABS(PA(2)-PB(2)).LT.EGEOM) THEN IF (ABS(PA(3)-PB(3)).LT.EGEOM) THEN WRITE(*,*) 'Vertices ',IV,JV,' are not distinct' LERROR=.TRUE. END IF END IF END IF 140 CONTINUE 130 CONTINUE IF (LERROR) THEN WRITE(*,*) WRITE(*,*) 'ERROR(LSEM3) - Vertices (see above) coincide' WRITE(*,*) 'Execution terminated' STOP END IF END IF C Check that the panels are not of disproportionate sizes IF (LVALID) THEN SIZMAX=ZERO SIZMIN=DIAM**2 DO 150 ISE=1,NPI QA(1)=VERTEX(PIELV(ISE,1),1) QA(2)=VERTEX(PIELV(ISE,1),2) QA(3)=VERTEX(PIELV(ISE,1),3) QB(1)=VERTEX(PIELV(ISE,2),1) QB(2)=VERTEX(PIELV(ISE,2),2) QB(3)=VERTEX(PIELV(ISE,2),3) QC(1)=VERTEX(PIELV(ISE,3),1) QC(2)=VERTEX(PIELV(ISE,3),2) QC(3)=VERTEX(PIELV(ISE,3),3) SIZE=AREA(QA,QB,QC) SIZMAX=MAX(SIZMAX,SIZE) SIZMIN=MIN(SIZMIN,SIZE) 150 CONTINUE IF (SIZMAX.GT.10.0D0*SIZMIN) THEN WRITE(*,*) 'WARNING(LSEM3) - panels of disproportionate' WRITE(*,*) ' sizes' END IF END IF C Validation of the surface functions IF (LVALID.AND.LSOL) THEN LERROR=.FALSE. DO 170 ISE=1,NPI IF (MAX(ABS(HA(ISE)),ABS(HB(ISE))).LT.1.0D-6) * LERROR=.TRUE. 170 CONTINUE IF (LERROR) THEN WRITE(*,*) WRITE(*,*) 'ERROR(LSEM3) - at most one of HA(i),HB(i)' WRITE(*,*) ' may be zero for all i' WRITE(*,*) 'Execution terminated' STOP END IF END IF C Compute the discrete L, M, Mt and N matrices C Loop(IPI) through the points on the shell DO 510 IPI=1,NPI C Set P PA(1)=VERTEX(PIELV(IPI,1),1) PA(2)=VERTEX(PIELV(IPI,1),2) PA(3)=VERTEX(PIELV(IPI,1),3) PB(1)=VERTEX(PIELV(IPI,2),1) PB(2)=VERTEX(PIELV(IPI,2),2) PB(3)=VERTEX(PIELV(IPI,2),3) PC(1)=VERTEX(PIELV(IPI,3),1) PC(2)=VERTEX(PIELV(IPI,3),2) PC(3)=VERTEX(PIELV(IPI,3),3) P(1)=(PA(1)+PB(1)+PC(1))/THREE P(2)=(PA(2)+PB(2)+PC(2))/THREE P(3)=(PA(3)+PB(3)+PC(3))/THREE C Set VECP to the normal on the shell of the panel at P CALL NORM3(PA,PB,PC,VECP) C Loop(IPI) through the panels DO 520 JPI=1,NPI C Set QA and QB, the coordinates of the edges of the ISEth panel QA(1)=VERTEX(PIELV(JPI,1),1) QA(2)=VERTEX(PIELV(JPI,1),2) QA(3)=VERTEX(PIELV(JPI,1),3) QB(1)=VERTEX(PIELV(JPI,2),1) QB(2)=VERTEX(PIELV(JPI,2),2) QB(3)=VERTEX(PIELV(JPI,2),3) QC(1)=VERTEX(PIELV(JPI,3),1) QC(2)=VERTEX(PIELV(JPI,3),2) QC(3)=VERTEX(PIELV(JPI,3),3) C Set LPONEL IF (IPI.EQ.JPI) THEN LPONEL=.TRUE. ELSE LPONEL=.FALSE. END IF C Select quadrature rule for L3LC C : Select the quadrature rule XQON-WQON in the case when the C : point p lies on the panel, otherwise select XQOFF-WQOFF C [Note that the overall method would benefit from selecting from C a wider set of quadrature rules, and an appropriate method C of selection] IF (LPONEL) THEN NQ=NQON DO 600 IQ=1,NQ XQ(IQ)=XQON(IQ) YQ(IQ)=YQON(IQ) WQ(IQ)=WQON(IQ) 600 CONTINUE ELSE NQ=NQOFF DO 610 IQ=1,NQ XQ(IQ)=XQOFF(IQ) YQ(IQ)=YQOFF(IQ) WQ(IQ)=WQOFF(IQ) 610 CONTINUE END IF C All operators are required LL=.TRUE. LM=.TRUE. LMT=.TRUE. LN=.TRUE. C Call L3LC CALL L3LC(P,VECP,QA,QB,QC,LPONEL, * MAXNQ,NQ,XQ,YQ,WQ, * LVAL,EGEOM,EQRULE,LFAIL1, * LL,LM,LMT,LN,DISL,DISM,DISMT,DISN) AMAT(IPI,JPI)=-DISM AMAT(IPI,NPI+JPI)=0.0D0 AMAT(NPI+IPI,JPI)=DISN AMAT(NPI+IPI,NPI+JPI)=0.0D0 BMAT(IPI,JPI)=-DISL BMAT(IPI,NPI+JPI)=0.0D0 BMAT(NPI+IPI,JPI)=DISMT BMAT(NPI+IPI,NPI+JPI)=0.0D0 C Close loop(JPI) 520 CONTINUE AMAT(IPI,NPI+IPI)=1.0D0 BMAT(NPI+IPI,NPI+IPI)=1.0D0 C Close loop(IPI) 510 CONTINUE C Set C,HALPHA,HBETA,F DO 700 IPI=1,NPI C Set C as the incident phi,v C(IPI)=HIPHI(IPI) C(NPI+IPI)=HIVEL(IPI) C Set HALPHA as the coefficient of [{\delta} {\Phi}] HALPHA(IPI)=HA(IPI) HALPHA(NPI+IPI)=HAA(IPI) C Set HBETA as the coefficient of [{\nu} V] HBETA(IPI)=HB(IPI) HBETA(NPI+IPI)=HBB(IPI) C Set F as [f F] HFFF(IPI)=HF(IPI) HFFF(NPI+IPI)=HFF(IPI) 700 CONTINUE C On solution C becomes {\zeta}/{\epsilon} at the central points CALL GLS(2*MAXNPI,2*NPI,AMAT,BMAT,C,HALPHA,HBETA,HFFF, * PHIINF,VELINF,LFAIL,PERM,XORY,WKSPC1,WKSPC2,WKSPC3) C Set C DO 710 IPI=1,NPI PHIDIF(IPI)=PHIINF(IPI) PHIAV(IPI)=PHIINF(NPI+IPI) VELDIF(IPI)=VELINF(IPI) VELAV(IPI)=VELINF(NPI+IPI) 710 CONTINUE C SOLUTION IN THE DOMAIN C Compute sound pressures at the selected exterior points. C Loop through the the points in the exterior region DO 800 IPE=1,NPE C Set P P(1)=PEXT(IPE,1) P(2)=PEXT(IPE,2) P(3)=PEXT(IPE,3) C Set VECP, this is arbitrary VECP(1)=ONE VECP(2)=ZERO VECP(3)=ZERO C Initialise SUMPHI to the incident potential SUMPHI=PEIPHI(IPE) C Loop(IPI) through the panels DO 850 JPI=1,NPI C Set QA and QB, the coordinates of the edges of the JPIth panel QA(1)=VERTEX(PIELV(JPI,1),1) QA(2)=VERTEX(PIELV(JPI,1),2) QA(3)=VERTEX(PIELV(JPI,1),3) QB(1)=VERTEX(PIELV(JPI,2),1) QB(2)=VERTEX(PIELV(JPI,2),2) QB(3)=VERTEX(PIELV(JPI,2),3) QC(1)=VERTEX(PIELV(JPI,3),1) QC(2)=VERTEX(PIELV(JPI,3),2) QC(3)=VERTEX(PIELV(JPI,3),3) C All the points do not lie on the shell hence LPONEL=.FALSE. LPONEL=.FALSE. C Only L,M operators are required. Set LL,LM true, LMT,LN false. LL=.TRUE. LM=.TRUE. LMT=.FALSE. LN=.FALSE. C Call L3LC. CALL L3LC(P,VECP,QA,QB,QC,LPONEL, * MAXNQ,NQ,XQ,YQ,WQ, * LVAL,EGEOM,EQRULE,LFAIL, * LL,LM,LMT,LN,DISL,DISM,DISMT,DISN) C Accumulate phi L_EH(IPE,IPI)=DISL M_EH(IPE,IPI)=DISM SUMPHI=SUMPHI+DISM*PHIDIF(JPI)-DISL*VELDIF(JPI) C Close loop (JPI) through the panels 850 CONTINUE PEPHI(IPE)=SUMPHI C Close loop(IPE) through the exterior points 800 CONTINUE END C ---------------------------------------------------------------------- C Subroutines required for L3LC (not in file L3LC.FOR) C Subroutine for returning the square root. REAL*8 FUNCTION FNSQRT(X) REAL*8 X FNSQRT=SQRT(X) END